function m1deltae(latt,lattm) result(m1de)
use para
implicit none
double precision::m1de
double precision,external::ran1,cauchys,prod
double precision::sx,sy,sz,squ,sx1,sy1,sz1,squ1,norm,newspin1,newspin2,newspin3  !,sdelta
double precision,dimension(1:3)::spinsum
integer,intent(in)::latt,lattm
integer::i,j,k,l1,m,n,temp11,temp22
i=latt

!the new spin configration is identified and deltaspin is calculated in subroutine deltaspin
!which is then used to calculate the energy change. 
!phi and theta are evenly distributed 
!if(sflag=='flat')then
!newtheta(latt,lattm)=int(ran2()*(thenum+1))
!newphi(latt,lattm)=int(ran2()*(phinum+1))

if(sadap/=1)then
!phi and theta are evenly distributed 
squ=1.0
do while(squ>0.5)
sx=ran1(seed1)-0.5
sy=ran1(seed1)-0.5
sz=ran1(seed1)-0.5
squ=sqrt(sx*sx+sy*sy+sz*sz)
end do
   newspin(1,latt,lattm)=sx/squ
   newspin(2,latt,lattm)=sy/squ
   newspin(3,latt,lattm)=sz/squ

   dspin(:,latt,lattm)=newspin(:,latt,lattm)-spin(:,latt,lattm)
   spinsum(:) = newspin(:,latt,lattm)+spin(:,latt,lattm) 

else 
!if 'adapt' is set to 1 then spin configration is updated adaptively 
squ1=2.0
do while(squ1>1.0)
sx1=2.*ran1(seed1)-1.0
sy1=2.*ran1(seed1)-1.0
sz1=2.*ran1(seed1)-1.0
squ1=sqrt(sx1*sx1+sy1*sy1+sz1*sz1)
end do
sx1=sx1*sdelta
sy1=sy1*sdelta
sz1=sz1*sdelta

newspin1=spin(1,latt,lattm)+sx1
newspin2=spin(2,latt,lattm)+sy1
newspin3=spin(3,latt,lattm)+sz1

norm=sqrt(newspin1*newspin1+newspin2*newspin2+newspin3*newspin3)

   newspin(1,latt,lattm)=newspin1/norm
   newspin(2,latt,lattm)=newspin2/norm
   newspin(3,latt,lattm)=newspin3/norm



 dspin(:,latt,lattm)=newspin(:,latt,lattm)-spin(:,latt,lattm)

 spinsum(:) = newspin(:,latt,lattm)+spin(:,latt,lattm)
end if

m1de=0.0
do j=1,natom, 1
  do l1=-mmldx, mmldx,1
    do m=-mmldy,mmldy,1
      do n=-mmldz,mmldz,1
        if(n==0 .and. &
          &m==0 .and. &
          &l1==0 .and. &
          &j==lattm) then
        cycle
        else
        m1de=m1de+&
        &jexc(lattm,j,n,m,l1)*(&
        &prod(dspin(1,latt,lattm),&
        &spin(:,mld(n,m,l1,latt),j)))! +&
        end if
      end do ! n
    end do ! m
  end do ! l
end do !j

if(anis_flag==1)then
  m1de=m1de-prod(D4,spinsum)*prod(D4,dspin(1,latt,lattm))
end if

if(hflag==1)then  
   m1de=m1de+prod(dspin(1,latt,lattm),hfield)*magmom
end if


end function
